Cite this Article
Gao Wei, Yu Cheng, Yao Feng. Droplets breakup via a splitting microchannel. Chinese Physics B, 2020, 29(5): 054702  
Permissions
Droplets breakup via a splitting microchannel
Gao Wei1, Yu Cheng2, †, Yao Feng3
School of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138, USA
Department of Mechanical Engineering, University of Hawaii at Manoa, Honolulu, HI 96822, USA
Jiangsu Key Laboratory of Micro and Nano Heat Fluid Flow Technology and Energy Application, School of Environmental Science and Engineering, Suzhou University of Science and Technology, Suzhou 215009, China

 

† Corresponding author. E-mail: chengy8@hawaii.edu

Abstract

On the basis of a volume of fluid (VOF) liquid/liquid interface tracking method, we apply a two-dimensional model to investigate the dynamic behaviors of droplet breakup through a splitting microchannel. The feasibility and applicability of the theoretical model are experimentally validated. Four flow regimes are observed in the splitting microchannel, that is, breakup with permanent obstruction, breakup with temporary obstruction, breakup with tunnels, and non-breakup. The results indicate that the increase of the capillary number Ca provides considerable upstream pressure to accelerate the droplet deformation, which is favorable for the droplet breakup. The decrease of the droplet size contributes to its shape changing from the plug to the sphere, which results in weakening droplet deformation ability and generating the non-breakup flow regime.

PACS: ;47.55.dr;;82.70.Kj;;47.55.db;;68.05.-n;
Keyword:;microfluidics;;droplet splitting;;hydrodynamic;;numerical simulation;
1. Introduction

Droplet-based microfluidics has attracted much attention in various fields, including chemical engineering,[1–3] bioengineering,[4,5] material science,[6] energy,[7–9] and beyond,[10,11] owing to its distinguish capacity in generating monodisperse droplets. The microchannels have significant flow resistance, which limits the fluid injecting rate in the range of several milliliters per hour.[12,13] It has to take prolonged time to achieve massive production of droplets, and remains a great challenge to enlarge the applications of droplets in the real environment.[14–16] Therefore, it is of particular importance to increase the droplet production rate and fully understand the hydrodynamic mechanisms to realize the controllable adjustment of the droplet behaviors.[17–19]

As an alternative method, the droplet splitting method has been proposed to increase the droplet production rate. Several theoretical and experimental studies focusing on the hydrodynamics process in different splitting microchannels have been carried out. Fu et al.[20] studied the bubble breakup process in the microfluidic T-junction considering the effect of the hydrodynamic feedback at the downstream microchannels. Carlson et al.[21] established a three-dimensional model based on the phase-field theory to study the droplet splitting process in a Y-junction numerically. It was demonstrated that the initial droplet size and capillary number Ca have significant influences on the flow regime, and the non-breakup flow regime causes the severe flow asymmetry once the droplet becomes close to the obstruction (Ca = μ V/σ, where μ and V are the viscosity and the velocity of the continuous phase, respectively and σ is the interfacial tension coefficient between the two fluid phases). Liu et al.[22] have developed a multiphase lattice Boltzmann model to simulate the contact-line dynamics, and applied this method to simulate the droplet breakup process in a microfluidic T-junction. They found that the droplet can break up into two unequal-sized daughter droplets in the non-ideal branch. Furthermore, Lee et al.[23] have demonstrated the droplets splitting behaviors in obstacles via numerical methods. It was found that the droplet has three modes, that is, the breakup, non-breakup, and breakup first and then merging. Besides, the use of long splitting obstacles can promote uniform daughter droplets after splitting.

Apart from the theoretical approaches, experimental studies have also been carried out. Liu and Zhang[24] developed a two-dimensional (2D) lattice Boltzmann model and studied the droplet formation in the T-junction. The effects of the capillary number, flow rate ratio, viscosity ratio, and contact angles were studied systematically. Chen and Deng[25] developed a phase-field multiphase lattice Boltzmann model and numerically investigated the hydrodynamic behaviors of a droplet passing through a microfluidic T-junction. The vortex flow inside droplets determined the breakup regime of the droplet. Link et al.[26] employed a three-level T-junction array to drive droplet split into three times in a row. Although this T-junction array could efficiently achieve droplet splitting, the space utilization was limited. While Abate and Weitz[27] utilized a parallel array to split large droplets into 1/16 of their original size, and the maximum production could reach 7000 μL/h. Moreover, Chen et al.[28] have developed glass capillary microfluidic devices to split both single emulsions and double emulsions, which successfully generated highly monodisperse droplets. However, the critical splitting mechanisms in this splitting microfluidic device are indistinct.

Herein, a numerical model based on the volume of fluid (VOF) liquid/liquid interface tracking method is applied to understand the droplet behaviors in a splitting microchannel. The flow regimes of droplets, including breakup with permanent obstruction, breakup with temporary obstruction, breakup with tunnels, and non-breakup, are studied with detailed pressure and velocity distributions. The effects of the droplet size and capillary number Ca are also discussed. Furthermore, a phase diagram clarifies the boundary between four flow regimes.

2. Theoretical model
2.1. Physical model

A 2D model is applied to investigate the droplets breakup behaviors in a splitting microchannel (see Figs. 1(a) and 1(b)). As illustrated in Fig. 1(c), the droplets and continuous phase are injected into the main channel, and the two branches of the splitting channel, that is, channels i and ii, are the same. The lengths of the main and splitting channels are lm = 2000 μm and ls = 3000 μm, respectively. The widths of the main and splitting channels are w = 300 μm and ws = 180 μm, respectively. The width of the splitting obstruction is wo = 60 μm. Besides, l represents the length of the droplet, and δ represents the neck of the splitting droplet. The continuous phase in the splitting microchannel is 2 wt% polyvinyl alcohol aqueous solution (PVA), and the droplet phase is ethoxylated trimethylolpropane triacrylate (ETPTA) oil. The PVA aqueous solution has a density of 1020 kg/m3 and a viscosity of 3 mPa⋅s. The ETPTA oil has a density of 1011 kg/m3 and a viscosity of 65 mPa⋅s. The interfacial tension between the droplet and continuous phase is 2.24 mN/m, and the wall of the capillary is hydrophilic with a contact angle of 70°.

Fig. 1. The geometric structure and physical model of the splitting microchannel: (a) microfluidic chip, (b) optical image of the splitting microchannel, (c) physical model.
2.2. Mathematical model

The splitting behavior of the droplet is a visible result of the liquid/liquid phase interface motion and deformation. Therefore, based on the physical model of droplet splitting, a two-dimensional splitting microchannel model is established, and a VOF method[29] is utilized to track the evolution and development of the interfaces.

2.2.1. Governing equation

In the VOF method, a volume fraction function αd is utilized to represent the proportion of droplet in one simulation cell as[30]

The sum of the volume fractions of the droplet phase, αd, and continuous phase, αc, in a cell is 1,

The conservation equations of the volume fraction can be express as

where t is the time, and U is the velocity. The U is governed by the continuum equation and Navier–Stokes equations and can be written as

where P is the pressure, ρ is the density, μ is the viscosity, and Fv is the source term. They can express as

where ρd and ρc are the densities of the droplet and continuous phase, μd and μc are the viscosities of the droplet and continuous phase, g is the gravity acceleration, and Fs is the volumetric surface tension force.

In the microchannel, the effect of gravity on the droplet behavior is much less than that of the interfacial tension; thus g can be neglected. The surface tension Fs is described by the continuous surface force (CSF) model. In the CSF model, Fs is thought as a bulk force acting on a cell of the phase interface region, which can be characterized as a function of interfacial tension coefficient and interface shape[31]

where σ is the interfacial tension coefficient, and κ is the mean curvature of the interface.

2.2.2. Boundary condition

As illustrated in Fig. 1(c), the origin of the coordinates is at the center of the inlet of the main channel. On the left of the computational domain, the boundary condition is the velocity entry

On the right of the computational domain, the boundary condition is the pressure outlet. The droplets are injected into the main channel until the continuous phase reaches the steady-state in the computational domain. The walls of the computational domain are no slip wall boundary. The VOF method includes an implicit slip length at no-slip boundaries which can explicitly track an interface that requires a slip condition for the mesh node at the contact line.[32,33]

2.2.3. Numerical methods

The finite volume method, combining with the finite difference method, is used to solve the governing equations numerically. The laminar model is adopted since the Reynold number is very low in the microchannel. The physical parameters of the fluid, such as density, viscosity, and interfacial tension coefficient, are constants. The liquid–liquid interface is reconstructed by the Geo–Reconstruct format of the sectional interface (PLIC) method. The semi-implicit method for the pressure linked equations (SIMPLE) algorithm is used for the pressure–velocity coupling interpolation. The body force weighted scheme is chosen for the pressure interpolation. The momentum equation is discretized by using a second-order upwind scheme. To achieve quick convergence, the under-relaxation factors implemented are 0.3 (pressure), 0.5 (density), 0.8 (source term), and 0.7 (momentum). The iteration in a one-time step ends when the relative residuals are less than 1%.

A mesh independence study is conducted using four different mesh numbers, that is, 306092, 349680, 361000, and 473570, respectively. As illustrated in Fig. 2, the phase interface morphologies of droplets splitting at the four different grid numbers are compared. When the mesh number is 349680, the degree of deformation of the droplets is very similar to the numerical results of 361000 and 473570. Also, the mesh density near the splitting obstruction is sufficient to track the deformation behaviors of the phase interface. Therefore, the total computational units, 349680, are used in the calculation.

Fig. 2. Grid independence test.
2.3. Model validation

In order to verify the feasibility of the established theoretical model, a series of verification experiments based on a high-speed microscopic visualization system are carried out. As shown in Fig. 3, the high-speed microscopic visualization system is mainly composed of the syringe pump unit, microfluidic chip, and high-speed microscopic imaging unit. The syringe pump unit injects the dispersed and continuous phase fluid into the microfluidic chip. The high-speed microscopic imaging unit can capture and record the droplet evolution processes in real-time.

Fig. 3. Experiment setup.

The numerical simulation results are utilized to compare with the experimental data to verify the rationality of the proposed theoretical model. In this case, the dimensionless parameters of time and neck thickness are used. The dimensionless time is tv = (t – t0)/(t1 – t0), where t0 is the initial time, t is the current time, and t1 is the end time. The dimensionless neck thickness is δ* = δ/w, where δ is the neck thickness, and w is the width of the main channel. The comparison among the three-dimensional (3D) experiment, the 2D experiment, and the 2D numerical simulations of droplet interface evolution in the splitting microchannel is performed in Fig. 4. As shown in Fig. 4, the 2D experiment results of droplet interface evolution and the dimensionless neck thickness evolution are similar to the 3D experiment results within the acceptable error. Moreover, as shown in Fig. 4(d), the 2D numerical results are in good agreement with the experimental results, except for a few differences at the end of the droplet breakup. Similar interfaces prove the correctness of the mathematical model. Thus, the 2D numerical simulation is adopted in the current study considering the computational efficiency and reasonability.

Fig. 4. The comparison of real-time image and numerical calculation results (Ca = 0.06, l/w = 2.35): (a) 3D real-time image of droplet splitting; (b) 2D real-time image of droplet splitting; (c) numerical results; (d) correlation of dimensionless neck thickness and dimensionless time.
3. Results and discussion
3.1. Flow regimes

Four typical flow regimes in the splitting microchannels are successfully reconstructed, that is, breakup with permanent obstruction, breakup with temporary obstruction, breakup with tunnels, and non-breakup, as demonstrated in Fig. 5. In the flow regime I, the droplet will block the main channel and the splitting channel all the time (see Fig. 5(a)). In the flow regime II, the droplets only block the splitting slit and do not contact with the wall of the splitting channel (see Fig. 5(b)). In the flow regime III, there is always a tunnel between the droplet and the wall to facilitate the continuous phase fluid (see Fig. 5(c)). The droplets breakup will not appear in the flow regime IV (see Fig. 5(d)). This geometrical division only considers the relationship between the droplets and the channel in the flow direction, and ignores the cross-sectional information of the droplets.

Fig. 5. Four typical flow regimes (Ca = 0.004): (a) breakup with permanent obstruction (l/w = 1.26); (b) breakup with temporary obstruction (l/w = 0.87); (c) breakup with tunnels (l/w = 0.70); (d) non-breakup (l/w = 0.62).
3.1.1. The breakup with permanent obstruction flow regime

The breakup with permanent obstruction flow regime can be divided into three stages, including entering, squeezing, and post-breakup, as illustrated in Fig. 6(a). The droplet first undergoes the entering stage in the main channel. At the initial time, the droplet blocks the main channel and gradually approaches to the splitting obstruction. After contacting the obstruction, its head enters the two branches of the splitting channel symmetrically and starts the squeezing stage. At the initiative of the squeezing stage, only the droplet head blocks the slit with a cashew nut shape, and most of the droplet still blocks the main channel (t* = 1.92 ms). Under the action of the interfacial tension, the rear of the droplet keeps convex in the main channel until most of the droplet enters the splitting microchannel (t* = 1.92–2.4 ms). At the end of the squeezing stage, the droplet head grows to block the splitting microchannel, and the droplet neck becomes thin near the splitting obstruction (t* = 2.8 ms). Finally, the droplet breaks into two individual daughter droplets (t* = 3.28–3.68 ms).

Fig. 6. Droplet behaviors for the breakup with permanent obstruction (Ca = 0.004, l/w = 1.26). Evolution of (a) droplet breakup morphology, (b) pressure field, and (c) velocity field.

As shown in Fig. 7(a), Poutlet i and Poutlet ii almost keep constant in the whole process. In the entering stage, the droplet flows from the main channel to the obstruction under the action of the upstream pressure, whose front and rear interfaces are undeformed. Thus, the upstream pressure Pinlet does not change in this stage. In the post-breakup stage, the droplet breaks into two small daughter droplets, resulting in a decrease of the flow resistance, which leads to the decrease of Pinlet.

Fig. 7. Evolution of the pressure for the breakup with permanent obstruction (Ca = 0.004, l/w = 1.26). (a) The inlet pressure of the main channel, Pinlet, and the outlet pressures of the splitting channel, Poutlet i and Poutlet ii. (b) The neck thickness δ*, the Laplace pressure of rear droplet interface Pδ, and the inlet and the outlet pressures in the A-A stage (t* = 0–3.60 ms).

The change in upstream pressure Pinlet in the squeezing stage is more complicated than that in other stages. As shown in Fig. 7(b), in the initial of the squeezing stage, the curvature radius of the rear increases slightly, while the front interface undergoes massive deformation and the curvature radius decreases, which promotes the pressure difference between the front and the rear, resulting in the increase of Pinlet until the droplet head blocks the slit. In the middle of squeezing, the pressure difference between the front and the rear decreases first and then increases, thus Pinlet decreases first and then increases. At the end of the squeezing stage, the curvature radius of the rear decreases, but the curvature radius of the front is almost constant, which cause a decrease in the pressure difference between the front and the rear, resulting in the decrease in the upstream pressure Pinlet. When Pinlet reaches the minimum, the curvature radius of the rear is almost constant, while the front continues to stretch under the action of the upstream pressure, which cause the pressure difference to increase, resulting in the slight increase in Pinlet.

3.1.2. The breakup with temporary obstruction flow regime

The breakup with temporary obstruction flow regime also has the same three stages as the flow regime I, as illustrated in Fig. 8. In the squeezing stage, the droplet head firstly enters two branches of the splitting channel symmetrically and blocks the slit with cashew nut shape (t* = 1.52 ms). Subsequently, the droplet still blocks the slit while its head has no contact with the wall until most of the droplet enters the splitting channel (t* = 1.52–1.787 ms). At the end of the squeezing stage, the slit is open, and the droplet head keeps away from the wall (t* = 2.03 ms). Finally, the droplet breaks up into two individual daughter droplets (t* = 2.48–2.88 ms).

Fig. 8. Droplet behaviors for the break up with temporary obstruction (Ca = 0.004, l/w = 0.87). Evolution of (a) droplet breakup morphology, (b) pressure field, and (c) velocity field.

The evolution of Poutlet i, Poutlet ii, and Pinlet in the entering stage in flow regime II is the same as that in the flow regime I (see Fig. 9). In the initial squeezing stage, Pinlet increases until the droplet head blocks the slit. In the middle of the squeezing stage, Pinlet decreases first and then increases. However, due to the curvature radius of the droplet does not change significantly, the fluctuation of Pinlet is small. At the end of the squeezing stage, the evolution of upstream pressure Pinlet is the same as that in the flow regime I, which shows a trend of decreasing first and then increasing.

Fig. 9. Evolution of the pressure for the breakup with temporary obstruction (Ca = 0.004, l/w = 0.87). (a) The inlet pressure of the main channel, Pinlet, and the outlet pressures of the splitting channel, Poutlet i and Poutlet ii. (b) The neck thickness δ*, the Laplace pressure of rear droplet interface Pδ, and the inlet and the outlet pressures in A-A stage (t* = 0–2.56 ms).
3.1.3. The breakup with tunnels flow regime

The breakup with tunnels flow regime also has the same three stages as the flow regimes I and II, as illustrated in Fig. 10. At the initial time, the spherical droplet flows and approaches to the obstruction without any contact with the wall. Subsequently, its head enters two branches of the splitting channel symmetrically with cashew nut shape (t* = 1.12 ms). In the whole squeezing stage, the droplets have no contact with the slit and wall all the time (t* = 1.12–1.36 ms). Finally, the droplet breaks into two individual daughter droplets (t* = 1.36–3.68 ms).

Fig. 10. Droplet behaviors for the breakup with tunnels (Ca = 0.004, l/w = 0.70). Evolution of (a) droplet breakup morphology, (b) pressure field, (c) velocity field.

The evolution of Poutlet i, Poutlet ii, Pinlet in the entering stage and post-breakup stage in flow regime III is the same as that in the flow regimes I and II (see Fig. 11). Since the droplet is small and has no contact with the wall in flow regime III, the dimensional limitation of the droplet is low, and the tunnel effect facilitates the continuous phase fluid flow conspicuously in the tunnel. Therefore, in the initial of the squeezing stage, when the droplet is in contact with the obstruction, the curvature radius of the rear increases which causes the decrease of the Laplace pressure at the rear interface, and the front interface undergoes deformation, resulting in the curvature radius decreasing, which contributes to increase the Laplace pressure at the front. Thus, the pressure difference between the front and the rear increases, increasing the upstream pressure Pinlet. In the middle of the squeezing stage, the droplet has been completely deformed. The curvature radius of the rear is reduced, resulting in an increase of the Laplace pressure at the rear. However, the Laplace pressure at the front is less changed, so the pressure difference between the front and the rear interfaces decreases, which causes a decrease of upstream pressure Pinlet.

Fig. 11. Evolution of the pressure for the breakup with tunnels (Ca = 0.004, l/w = 0.70). (a) The inlet pressure of the main channel, Pinlet, and the outlet pressures of the splitting channel, Poutlet i and Poutlet ii. (b) The neck thickness δ*, the Laplace pressure of rear droplet interface Pδ, and the inlet and the outlet pressures in A-A stage (t* = 0–3.5 ms).
3.1.4. The non-breakup flow regime

The non-breakup flow regime can be divided into three stages, including entering, sliding, and non-breakup (see Fig. 12). The entering stage (t* = 0–0.267 ms) is similar to that in flow regime III. After the front of the droplet contacts with the splitting obstruction, the droplet is compressed into a cashew nut shape by the upstream pressure (t* = 1.09 ms). Subsequently, the front of the cashew nut shape will select one branch of the splitting channel for orientational motion (t* = 1.09–1.28 ms). Finally, the unbalanced droplet passes over the obstruction and quickly enters one branch without breakup.

Fig. 12. Droplet behaviors for the non-breakup flow (Ca = 0.004, l/w = 0.62). Evolution of (a) droplet non-breakup morphology, (b) pressure field, and (c) velocity field.

As illustrated in Fig. 13, in the entering stage, Pinlet, Poutlet i, and Poutlet ii are almost unchanged (t* = 0–0.267 ms). In the sliding stage, the inlet pressure Pinlet continues to increase and reaches the maximum at t* = 1.1 ms when the rear of the droplet is closest to the wall. Before t* = 1.1 ms, Poutlet i and Poutlet ii begin to become imbalance, which shows an increase of Poutlet i and a decrease of Poutlet ii. Meanwhile, the droplet begins to select the sliding direction. When Pinlet reaches the maximum, the pressure imbalance continues to expand, resulting in a significant bias in the direction of the droplet movement. Finally, the droplet loses balance and passes through the splitting obstruction to enter the branch via the slit. Moreover, a vortex that favors the deformation of the droplet appears at the rear of the droplet during the sliding stage.

Fig. 13. Evolution of the pressure for the non-breakup flow (Ca = 0.004, l/w = 0.62). (a) The inlet pressure of the main channel, Pinlet, and the outlet pressures of the splitting channel, Poutlet i and Poutlet ii. (b) The neck thickness δ*, the Laplace pressure of rear droplet interface Pδ, and the inlet and the outlet pressures in A-A stage (t* = 0–1.6 ms).
3.2. Phase diagram of droplet splitting

As seen from Fig. 14, when the dimensionless droplet length l/w is shorter than 0.8, the flow regime will change from non-breakup to breakup with tunnels as the capillary number Ca increases. Large capillary numbers can provide considerable higher upstream pressure and promote the deformation of the droplets so that the increase in the capillary number facilitates the occurrence of droplet breakup. When the dimensionless length l/w is greater than 1, high capillary numbers make the breakup process change from permanent obstruction to temporary obstruction due to the acceleration of droplets passing through the splitting channel. In addition, the flow regimes will change from breakup to non-breakup as the dimensionless length l/w decreases. When the dimensionless length reduces, the droplet shape changes from plug to the sphere, which causes contact with the wall surface to reduce. Thus, the transition rate of the droplet interface decreases and non-breakup eventually happens. Leshansky and Pismen[34] first proposed the equation to describe the boundary between the breakup and non-breakup. The correlation between the dimensionless droplet length and Ca is written as

where l* = l/w, and m and n are constant in this condition. After comparing the present numerical simulation results and fitting the curve with the boundary, we can conclude that the predict equation is

As shown in Fig. 14, the fitting curve represents the boundary between the breakup and non-breakup, which shows a satisfying agreement.

Fig. 14. Phase diagram of droplet flow in a splitting microchannel.
4. Conclusion and perspectives

Based on the VOF liquid/liquid interface tracking method, a two-dimensional model is applied to investigate the dynamic behaviors of droplet breakup via a splitting microchannel. The feasibility and applicability of the theoretical model are experimentally verified. Numerical simulation is carried out to study the droplet splitting process and reveal the evolution process of the pressure field and the velocity field during the droplet breakup. The breakup mechanisms of the droplet in different flow regimes are obtained. The conclusions can be summarized as follows.

(1) Four flow regimes are observed in the splitting microchannel, including breakup with permanent obstruction, breakup with temporary obstruction, breakup with tunnels, and non-breakup. Three stages of the breakup process are entering, squeezing, and post-breakup. The three stages of the non-breakup process are entering, sliding, and non-breakup.

(2) In the entering stage, the front and rear interfaces of the droplet are almost uniform under the action of the upstream pressure, which contributes to a constant inlet pressure Pinlet of the main channel. In the post-breakup stage, the droplet is split into two small portions in the splitting microchannel, resulting in the decrease of the flow resistance, which leads to a rapid decrease of Pinlet.

(3) In the squeezing stage, the droplet is driven by the upstream pressure and hindered by the splitting structure. During the deformation process, the Laplace pressure difference between the front and rear interfaces of the droplet is positively correlated with Pinlet.

(4) The increase of the capillary number Ca provides a more considerable upstream pressure to promote the deformation of the droplet, which is favorable for the breakup of the droplet. The decrease of the droplet size helps to change the shape of the droplet from plug to the sphere, resulting in the reduction of the droplet deformation and the appearance of the non-breakup flow regime.

Reference
[1] Zhang C Shen C Chen Y 2017 Int. J. Heat Mass Transfer 104 1135
[2] Vladisavljevic G T Khalid N Neves M A Kuroiwa T Nakajima M Uemura K Ichikawa S Kobayashi I 2013 Adv. Drug Deliv Rev. 65 1626
[3] Li J Liu H Ioannou N Zhang Y Reese J M 2015 Commun. Comput. Phys. 17 1113
[4] Shang L Cheng Y Zhao Y 2017 Chem. Rev. 117 7964
[5] Li Y Jain M Ma Y Nandakumar K 2015 Soft Matter 11 3884
[6] Fu F Chen Z Zhao Z Wang H Shang L Gu Z Zhao Y 2017 Proc. Natl Acad. Sci. USA 114 5900
[7] Lan K Liu J Li Z et al. 2016 Matter Radiat. Extremes 1 8
[8] Zhang C Yu F Li X Chen Y 2019 AlChE J. 65 1119
[9] Campbell E M Goncharov V N Sangster T C et al. 2017 Matter Radiat. Extremes 2 37
[10] Lee T Y Choi T M Shim T S Frijns R A Kim S H 2016 Lab Chip 16 3415
[11] Liang W G Yang C Wen G Q Wang W Ju X J Xie R Chu L Y 2014 Appl. Therm. Eng. 70 817
[12] Chen Y Zhang C Shi M Yang Y 2009 AlChE J. 56 2018
[13] Wang N Liu H Zhang C 2017 J. Rheol. 61 741
[14] Ofner A Mattich I Hagander M Dutto A Seybold H Rühs P A Studart A R 2019 Adv. Funct. Mater. 29 1806821
[15] Nisisako T Ando T Hatsuzawa T 2012 Lab Chip 12 3426
[16] Hoang D A Haringa C Portela L M Kreutzer M T Kleijn C R van Steijn V 2014 Chem. Eng. J. 236 545
[17] Chen Y Liu X Shi M 2013 Appl. Phys. Lett. 102 051609
[18] Zhang C Deng Z Chen Y 2014 Int. J. Heat Mass Transfer 70 322
[19] Liu H Ba Y Wu L Li Z Xi G Zhang Y 2018 J. Fluid Mech. 837 381
[20] Fu T Ma Y Li H Z 2014 AlChE J. 60 1920
[21] Carlson A Do-Quang M Amberg G 2010 Int. J. Multiphase Flow 36 397
[22] Liu H Ju Y Wang N Xi G Zhang Y 2015 Phys. Rev. 92 033306
[23] Lee J Lee W Son G 2013 J. Mech. Sci. Technol. 27 1693
[24] Liu H Zhang Y 2009 J. Appl. Phys. 106 034906
[25] Chen Y Deng Z 2017 J. Fluid Mech. 819 401
[26] Link D R Anna S L Weitz D A Stone H A 2004 Phys. Rev. Lett. 92 054503
[27] Abate A R Weitz D A 2011 Lab Chip 11 1911
[28] Chen Y Gao W Zhang C Zhao Y 2016 Lab Chip 16 1332
[29] Hirt C W Nichols B D 1981 J. Comput. Phys. 39 201
[30] Cristini V Tan Y C 2004 Lab Chip 4 257
[31] Brackbill J U Kothe D B Zemach C 1992 J. Comput. Phys. 100 335
[32] Afkhami S Zaleski S Bussmann M 2009 J. Comput. Phys. 228 5370
[33] Renardy M Renardy Y Li J 2001 J. Comput. Phys. 171 243
[34] Leshansky A M Pismen L M 2009 Phys. Fluids 21 023303